library(psych)
library(ggplot2)
library(tidyverse)
library(effects)
## Library for the multiple comparisons
library(phia)
## Library to compute the effect sizes
library(sjstats)
library(lsr)
library(lavaan)
library(ggthemes)
library(tibble)
library(car)
library(olsrr)
library(pwr)
library(R2jags)
library(rjags)
library(runjags)
library("lattice")
library("superdiag")
library(gtools)
library(bayestestR)
library("ggmcmc")
library(brms)
library(BayesFactor)
library(MCMCvis)
source('Useful Code V2.R')
library(readr)
attobo <- read_csv("attobo1_clean.csv")
colnames(attobo) <- tolower(colnames(attobo))
attobo$scenario <- as.factor(attobo$scenario)
attobo$attribution <- attobo$scenario
attobo$attribution_effort <- factor(attobo$attribution, levels = c("Effort", "Strategy", "Aptitude"))
attobo$attribution_strat <- factor(attobo$attribution, levels = c( "Strategy","Effort", "Aptitude"))
attobo$expectancy <- (attobo$expectancy - 1)/10
attobo$expectancy <- ifelse(attobo$expectancy == 1, .999,
ifelse(attobo$expectancy == 0, .001, attobo$expectancy) )
Descriptive Stats
## Sample Breakdown
group_by(attobo, attribution) %>%
summarise(n())
## means and sd of DVs
aggregate(cbind(combined_method, bruteforce, seekhelp)~attribution, data=attobo, FUN=mean)
aggregate(cbind(combined_method, bruteforce, seekhelp)~attribution, data=attobo, FUN=sd)
aggregate(cbind(total, convergent_method)~attribution, data=attobo, FUN=mean)
aggregate(cbind(total, convergent_method)~attribution, data=attobo, FUN=sd)
## Locus, Variability, and Uncontrollability
attobo %>%
group_by(attribution) %>%
summarise(meanLocusOut = mean(locusout, na.rm = T),
sdLocusOut = sd(locusout, na.rm = T),
meanVar = mean(variability, na.rm = T),
sdVar = sd(variability, na.rm = T),
meanUncon = mean(uncontrollability, na.rm = T),
sdUncon = sd(uncontrollability, na.rm = T))
attobo.advice <- attobo %>%
select(combined_method, bruteforce, seekhelp, convergent_method,
total, expectancy, total_method)
attobo.attribution <- attobo %>%
select(aptout, aptvar, aptuncon,
effortout, effortvar, effortuncon,
stratout, stratvar, stratuncon,
luckout, luckvar, luckuncon,
taskout, taskvar, taskuncon,
moodout, moodvar, mooduncon)
attobo.others <- attobo %>%
select(expectancy, gender, age, education, sms_mean,
smo_mean, mp_mean, so_mean)
ggplot(data = gather(attobo.advice, factor_key = TRUE), aes(x = factor(value))) +
geom_bar() +
facet_wrap(~ key,scales = "free", as.table = TRUE, nrow = 3) + xlab("") +
theme_bw()

ggplot(data = gather(attobo.attribution, factor_key = TRUE), aes(x = factor(value))) +
geom_bar() +
facet_wrap(~ key,scales = "free", as.table = TRUE, nrow = 3) + xlab("") +
theme_bw()

ggplot(data = gather(attobo.others, factor_key = TRUE), aes(x = factor(value))) +
geom_bar() +
facet_wrap(~ key,scales = "free", as.table = TRUE, nrow = 3) + xlab("") +
theme_bw()

Frequentist Estimate
library(glm2)
library(betareg)
method.freq <- glm2(combined_method ~ attribution_strat, data = attobo, family = poisson(link = "log"))
convergent.freq <- glm2(convergent_method ~ attribution_effort, data = attobo, family = poisson(link = "log"))
bruteforce.freq <- glm2(bruteforce ~ attribution_effort, data = attobo, family = poisson(link = "log"))
expectancy.freq <- betareg(expectancy ~ attribution, data = attobo, link = "logit" )
seekhelp.freq <- glm2(seekhelp ~ attribution, data = attobo, family = poisson(link = "log"))
total.freq <- glm2(total ~ attribution, data = attobo, family = poisson(link = "log"))
summary(glm2(combined_method ~ attribution_strat, data = attobo, family = poisson(link = "log")))
##
## Call:
## glm2(formula = combined_method ~ attribution_strat, family = poisson(link = "log"),
## data = attobo)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.74895 -1.18818 0.09237 0.36315 2.69974
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.42488 0.09806 4.333 1.47e-05 ***
## attribution_stratEffort -0.77319 0.17449 -4.431 9.38e-06 ***
## attribution_stratAptitude -0.51870 0.16127 -3.216 0.0013 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 209.38 on 202 degrees of freedom
## Residual deviance: 186.37 on 200 degrees of freedom
## AIC: 508.8
##
## Number of Fisher Scoring iterations: 5
summary(glm2(convergent_method ~ attribution_effort, data = attobo, family = poisson(link = "log")))
##
## Call:
## glm2(formula = convergent_method ~ attribution_effort, family = poisson(link = "log"),
## data = attobo)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.4038 -0.8639 -0.5688 0.4310 1.8607
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.01482 0.12217 -0.121 0.903
## attribution_effortStrategy -1.80680 0.32528 -5.555 2.78e-08 ***
## attribution_effortAptitude -0.97100 0.23435 -4.143 3.42e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 209.45 on 202 degrees of freedom
## Residual deviance: 161.03 on 200 degrees of freedom
## AIC: 339.06
##
## Number of Fisher Scoring iterations: 5
summary(glm2(bruteforce ~ attribution_effort, data = attobo, family = poisson(link = "log")))
##
## Call:
## glm2(formula = bruteforce ~ attribution_effort, family = poisson(link = "log"),
## data = attobo)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.8911 -0.6229 -0.6229 -0.3430 2.6321
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.9237 0.1925 -4.800 1.59e-06 ***
## attribution_effortStrategy -1.9095 0.5358 -3.564 0.000365 ***
## attribution_effortAptitude -0.7161 0.3376 -2.121 0.033905 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 167.14 on 202 degrees of freedom
## Residual deviance: 147.77 on 200 degrees of freedom
## AIC: 226.05
##
## Number of Fisher Scoring iterations: 6
summary(betareg(expectancy ~ attribution, data = attobo, link = "logit" ))
##
## Call:
## betareg(formula = expectancy ~ attribution, data = attobo, link = "logit")
##
## Standardized weighted residuals 2:
## Min 1Q Median 3Q Max
## -4.6329 -0.4291 -0.1732 0.1347 3.3119
##
## Coefficients (mean model with logit link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.740427 0.130347 5.680 1.34e-08 ***
## attributionEffort 0.440632 0.181451 2.428 0.0152 *
## attributionStrategy -0.002054 0.179944 -0.011 0.9909
##
## Phi coefficients (precision model with identity link):
## Estimate Std. Error z value Pr(>|z|)
## (phi) 2.4961 0.2267 11.01 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Type of estimator: ML (maximum likelihood)
## Log-likelihood: 64.7 on 4 Df
## Pseudo R-squared: 0.03119
## Number of iterations: 19 (BFGS) + 2 (Fisher scoring)
summary(glm2(total ~ attribution, data = attobo, family = poisson(link = "log")))
##
## Call:
## glm2(formula = total ~ attribution, family = poisson(link = "log"),
## data = attobo)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.22288 -0.75504 0.03143 0.51507 2.34979
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.78574 0.08248 9.527 <2e-16 ***
## attributionEffort 0.11872 0.11294 1.051 0.293
## attributionStrategy -0.11490 0.11967 -0.960 0.337
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 128.65 on 202 degrees of freedom
## Residual deviance: 124.56 on 200 degrees of freedom
## AIC: 650.56
##
## Number of Fisher Scoring iterations: 4
Anova(glm2(combined_method ~ attribution_strat, data = attobo, family = poisson(link = "log")))
Anova(glm2(convergent_method ~ attribution_effort, data = attobo, family = poisson(link = "log")))
Anova(glm2(bruteforce ~ attribution_effort, data = attobo, family = poisson(link = "log")))
Anova(betareg(expectancy ~ attribution, data = attobo, link = "logit" ))
Anova(glm2(total ~ attribution, data = attobo, family = poisson(link = "log")))
Bayesian Model - method Advice
strategy <- ifelse(attobo$attribution == "Strategy", 1,0)
effort <- ifelse(attobo$attribution == "Effort", 1, 0)
method <- attobo$combined_method
attobo.datjags <- list(strategy = strategy, effort = effort, method = method, N = nrow(attobo))
attobo.datjags
Model Specification
attobo.model <- function() {
for(i in 1:N){
# lambda is mean and variance of the poisson distribution
method[i] ~ dpois(lambda[i])
log(lambda[i]) <- b[1] + b[2] * effort[i] + b[3] * strategy[i]
}
#priors
b[1] ~ dnorm(0, 0.0025)
b[2] ~ dnorm(0, 0.0025)
b[3] ~ dnorm(0, 0.0025)
for (i in 1:N){
error[i] <- lambda[i] - method[i]
}
}
attobo.params <- c("b", "lambda", "error")
b.inits <- method.freq$coefficients
attobo.inits1 <- list( "b" = b.inits)
attobo.inits2 <- list( "b" = rep(0, 3) )
attobo.inits3 <- list( "b" = rep(1, 3) )
attobo.inits4 <- list( "b" = rep(-1, 3) )
attobo.inits5 <- list( "b" = rep(2, 3) )
attobo.inits <- list(attobo.inits1, attobo.inits2,attobo.inits3, attobo.inits4,attobo.inits5)
set.seed(1128)
attobo.fit <- jags(data = attobo.datjags, inits = attobo.inits, parameters.to.save = attobo.params,
n.chains = 5, n.iter = 100000, n.burnin = 10000, model.file = attobo.model)
attobo.fit.upd <- update(attobo.fit, n.iter =1000)
attobo.fit.upd <- autojags(attobo.fit)
attobo.fit.mcmc <- as.mcmc(attobo.fit)
save(attobo.fit.mcmc, file = "attobo.method.mcmc.RData")
Convergence Diagnostics
load("attobo.method.mcmc.RData")
MCMCtrace(attobo.fit.mcmc[,1:4], ind = TRUE, pdf = FALSE )


autocorr.plot(attobo.fit.mcmc[1][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[2][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[3][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[4][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[5][,1:4], ask = FALSE)

attobo.out <- as.data.frame(as.matrix(attobo.fit.mcmc))
## saving data
save(attobo.out, file = "attobo.method.RData")
load("attobo.method.RData")
pred.method <- attobo.out[, grep("lambda[", colnames(attobo.out), fixed = T)]
## estimated count for each participant
pred.method <- pred.method[, c(mixedsort(names(pred.method)))]
strategy.id <- as.numeric(unlist(which(attobo$attribution == "Strategy")))
effort.id <- as.numeric(unlist(which(attobo$attribution == "Effort")))
aptitude.id <- as.numeric(unlist(which(attobo$attribution == "Aptitude")))
pred.method$strategy.sim <- apply(pred.method[, strategy.id], 1, mean)
pred.method$effort.sim <- apply(pred.method[, effort.id], 1, mean)
pred.method$aptitude.sim <- apply(pred.method[, aptitude.id], 1, mean)
quantile(pred.method$strategy.sim, probs = c(.025, .5, .975))
## 2.5% 50% 97.5%
## 1.250904 1.526093 1.835312
quantile(pred.method$effort.sim, probs = c(.025, .5, .975))
## 2.5% 50% 97.5%
## 0.5170229 0.7007571 0.9303260
quantile(pred.method$aptitude.sim, probs = c(.025, .5, .975))
## 2.5% 50% 97.5%
## 0.6974865 0.9049689 1.1477881
pred.method$diff.1 <- pred.method$strategy.sim - pred.method$aptitude.sim
pred.method$diff.2 <- pred.method$strategy.sim - pred.method$effort.sim
pred.method$diff.3 <- pred.method$effort.sim - pred.method$aptitude.sim
quantile(pred.method$diff.1, probs = c(.025, .5, .975))
## 2.5% 50% 97.5%
## 0.2514705 0.6191405 0.9899373
quantile(pred.method$diff.2, probs = c(.025, .5, .975))
## 2.5% 50% 97.5%
## 0.4776900 0.8261376 1.1801479
quantile(pred.method$diff.3, probs = c(.025, .5, .975))
## 2.5% 50% 97.5%
## -0.5085726 -0.2057398 0.1046953
error.method <- attobo.out[, grep("error[", colnames(attobo.out), fixed = T)]
error <- as.vector(unlist(lapply(error.method, mean)))
mean(error)
## [1] 0.001381898

Expected Counts for each attribution
attobo.posterior.coefs <- select(pred.method, `strategy.sim`, `effort.sim`,`aptitude.sim`)
names(attobo.posterior.coefs) <- c("Strategy", "Effort", "Aptitude")
attobo.posterior.coefs.long <- gather(attobo.posterior.coefs)
#head(attobo.posterior.coefs.long)
library("ggridges")
ggplot(data = attobo.posterior.coefs.long, aes(x = value, y = key)) +
stat_density_ridges(quantile_lines = TRUE,
quantiles = c(0.025, 0.5, 0.975),
alpha = 0.7)

Bayesian Model - bruteforce Advice
strategy <- ifelse(attobo$attribution == "Strategy", 1,0)
effort <- ifelse(attobo$attribution == "Effort", 1, 0)
bruteforce <- attobo$bruteforce
attobo.datjags <- list(strategy = strategy, effort = effort, bruteforce = bruteforce, N = nrow(attobo))
attobo.datjags
Model Specification
attobo.model <- function() {
for(i in 1:N){
# lambda is mean and variance of the poisson distribution
bruteforce[i] ~ dpois(lambda[i])
log(lambda[i]) <- b[1] + b[2] * effort[i] + b[3] * strategy[i]
}
#priors
b[1] ~ dnorm(0, 0.0025)
b[2] ~ dnorm(0, 0.0025)
b[3] ~ dnorm(0, 0.0025)
for (i in 1:N){
error[i] <- lambda[i] - bruteforce[i]
}
}
attobo.params <- c("b", "lambda", "error")
b.inits <- bruteforce.freq$coefficients
attobo.inits1 <- list( "b" = b.inits)
attobo.inits2 <- list( "b" = rep(0, 3) )
attobo.inits3 <- list( "b" = rep(1, 3) )
attobo.inits4 <- list( "b" = rep(-1, 3) )
attobo.inits5 <- list( "b" = rep(2, 3) )
attobo.inits <- list(attobo.inits1, attobo.inits2,attobo.inits3, attobo.inits4,attobo.inits5)
set.seed(1128)
attobo.fit <- jags(data = attobo.datjags, inits = attobo.inits, parameters.to.save = attobo.params,
n.chains = 5, n.iter = 100000, n.burnin = 10000, model.file = attobo.model)
attobo.fit.upd <- update(attobo.fit, n.iter =1000)
attobo.fit.upd <- autojags(attobo.fit)
attobo.fit.mcmc <- as.mcmc(attobo.fit)
save(attobo.fit.mcmc, file = "attobo.bruteforce.mcmc.RData")
Convergence Diagnostics
load("attobo.bruteforce.mcmc.RData")
MCMCtrace(attobo.fit.mcmc[,1:4], ind = TRUE, pdf = FALSE )


autocorr.plot(attobo.fit.mcmc[1][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[2][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[3][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[4][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[5][,1:4], ask = FALSE)

attobo.out <- as.data.frame(as.matrix(attobo.fit.mcmc))
## saving data
save(attobo.out, file = "attobo.bruteforce.RData")
pred.bruteforce <- attobo.out[, grep("lambda[", colnames(attobo.out), fixed = T)]
## estimated count for each participant
pred.bruteforce <- pred.bruteforce[, c(mixedsort(names(pred.bruteforce)))]
strategy.id <- as.numeric(unlist(which(attobo$attribution == "Strategy")))
effort.id <- as.numeric(unlist(which(attobo$attribution == "Effort")))
aptitude.id <- as.numeric(unlist(which(attobo$attribution == "Aptitude")))
pred.bruteforce$strategy.sim <- apply(pred.bruteforce[, strategy.id], 1, mean)
pred.bruteforce$effort.sim <- apply(pred.bruteforce[, effort.id], 1, mean)
pred.bruteforce$aptitude.sim <- apply(pred.bruteforce[, aptitude.id], 1, mean)
quantile(pred.bruteforce$strategy.sim, probs = c(.025, .5, .975))
## 2.5% 50% 97.5%
## 1.250904 1.526093 1.835312
quantile(pred.bruteforce$effort.sim, probs = c(.025, .5, .975))
## 2.5% 50% 97.5%
## 0.5170229 0.7007571 0.9303260
quantile(pred.bruteforce$aptitude.sim, probs = c(.025, .5, .975))
## 2.5% 50% 97.5%
## 0.6974865 0.9049689 1.1477881
error.bruteforce <- attobo.out[, grep("error[", colnames(attobo.out), fixed = T)]
error <- as.vector(unlist(lapply(error.bruteforce, mean)))
mean(error)
## [1] 0.001381898

Expected Counts for each attribution
attobo.posterior.coefs <- select(pred.bruteforce, `strategy.sim`, `effort.sim`,`aptitude.sim`)
names(attobo.posterior.coefs) <- c("Strategy", "Effort", "Aptitude")
attobo.posterior.coefs.long <- gather(attobo.posterior.coefs)
#head(attobo.posterior.coefs.long)
library("ggridges")
ggplot(data = attobo.posterior.coefs.long, aes(x = value, y = key)) +
stat_density_ridges(quantile_lines = TRUE,
quantiles = c(0.025, 0.5, 0.975),
alpha = 0.7)

Bayesian Model - mandatory Advice
strategy <- ifelse(attobo$attribution == "Strategy", 1,0)
effort <- ifelse(attobo$attribution == "Effort", 1, 0)
convergent <- attobo$convergent_method
attobo.datjags <- list(strategy = strategy, effort = effort, convergent = convergent, N = nrow(attobo))
attobo.datjags
Model Specification
attobo.model <- function() {
for(i in 1:N){
# lambda is mean and variance of the poisson distribution
convergent[i] ~ dpois(lambda[i])
log(lambda[i]) <- b[1] + b[2] * effort[i] + b[3] * strategy[i]
}
#priors
b[1] ~ dnorm(0, 0.0025)
b[2] ~ dnorm(0, 0.0025)
b[3] ~ dnorm(0, 0.0025)
for (i in 1:N){
error[i] <- lambda[i] - convergent[i]
}
}
attobo.params <- c("b", "lambda", "error")
b.inits <- convergent.freq$coefficients
attobo.inits1 <- list( "b" = b.inits)
attobo.inits2 <- list( "b" = rep(0, 3) )
attobo.inits3 <- list( "b" = rep(1, 3) )
attobo.inits4 <- list( "b" = rep(-1, 3) )
attobo.inits5 <- list( "b" = rep(2, 3) )
attobo.inits <- list(attobo.inits1, attobo.inits2,attobo.inits3, attobo.inits4,attobo.inits5)
set.seed(1128)
attobo.fit <- jags(data = attobo.datjags, inits = attobo.inits, parameters.to.save = attobo.params,
n.chains = 5, n.iter = 100000, n.burnin = 10000, model.file = attobo.model)
attobo.fit.upd <- update(attobo.fit, n.iter =1000)
attobo.fit.upd <- autojags(attobo.fit)
attobo.fit.mcmc <- as.mcmc(attobo.fit)
save(attobo.fit.mcmc, file = "attobo.convergent.mcmc.RData")
Convergence Diagnostics
load("attobo.convergent.mcmc.RData")
MCMCtrace(attobo.fit.mcmc[,1:4], ind = TRUE, pdf = FALSE )


autocorr.plot(attobo.fit.mcmc[1][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[2][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[3][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[4][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[5][,1:4], ask = FALSE)

attobo.out <- as.data.frame(as.matrix(attobo.fit.mcmc))
## saving data
save(attobo.out, file = "attobo.convergent.RData")
load("attobo.convergent.RData")
pred.convergent <- attobo.out[, grep("lambda[", colnames(attobo.out), fixed = T)]
## estimated count for each participant
pred.convergent <- pred.convergent[, c(mixedsort(names(pred.convergent)))]
strategy.id <- as.numeric(unlist(which(attobo$attribution == "Strategy")))
effort.id <- as.numeric(unlist(which(attobo$attribution == "Effort")))
aptitude.id <- as.numeric(unlist(which(attobo$attribution == "Aptitude")))
pred.convergent$strategy.sim <- apply(pred.convergent[, strategy.id], 1, mean)
pred.convergent$effort.sim <- apply(pred.convergent[, effort.id], 1, mean)
pred.convergent$aptitude.sim <- apply(pred.convergent[, aptitude.id], 1, mean)
quantile(pred.convergent$strategy.sim, probs = c(.025, .5, .975))
## 2.5% 50% 97.5%
## 0.08305286 0.15728183 0.26887682
quantile(pred.convergent$effort.sim, probs = c(.025, .5, .975))
## 2.5% 50% 97.5%
## 0.7675159 0.9802979 1.2365364
quantile(pred.convergent$aptitude.sim, probs = c(.025, .5, .975))
## 2.5% 50% 97.5%
## 0.2390818 0.3637082 0.5336015
error.convergent <- attobo.out[, grep("error[", colnames(attobo.out), fixed = T)]
error <- as.vector(unlist(lapply(error.convergent, mean)))
mean(error)
## [1] 0.000883394

Expected Counts for each attribution
attobo.posterior.coefs <- select(pred.convergent, `strategy.sim`, `effort.sim`,`aptitude.sim`)
names(attobo.posterior.coefs) <- c("Strategy", "Effort", "Aptitude")
attobo.posterior.coefs.long <- gather(attobo.posterior.coefs)
#head(attobo.posterior.coefs.long)
library("ggridges")
ggplot(data = attobo.posterior.coefs.long, aes(x = value, y = key)) +
stat_density_ridges(quantile_lines = TRUE,
quantiles = c(0.025, 0.5, 0.975),
alpha = 0.7)

Bayesian Model - seekhelp Advice
strategy <- ifelse(attobo$attribution == "Strategy", 1,0)
effort <- ifelse(attobo$attribution == "Effort", 1, 0)
seekhelp <- attobo$seekhelp
attobo.datjags <- list(strategy = strategy, effort = effort, seekhelp = seekhelp, N = nrow(attobo))
attobo.datjags
Model Specification
attobo.model <- function() {
for(i in 1:N){
# lambda is mean and variance of the poisson distribution
seekhelp[i] ~ dpois(lambda[i])
log(lambda[i]) <- b[1] + b[2] * effort[i] + b[3] * strategy[i]
}
#priors
b[1] ~ dnorm(0, 0.0025)
b[2] ~ dnorm(0, 0.0025)
b[3] ~ dnorm(0, 0.0025)
for (i in 1:N){
error[i] <- lambda[i] - seekhelp[i]
}
}
attobo.params <- c("b", "lambda", "error")
b.inits <- seekhelp.freq$coefficients
attobo.inits1 <- list( "b" = b.inits)
attobo.inits2 <- list( "b" = rep(0, 3) )
attobo.inits3 <- list( "b" = rep(1, 3) )
attobo.inits4 <- list( "b" = rep(-1, 3) )
attobo.inits5 <- list( "b" = rep(2, 3) )
attobo.inits <- list(attobo.inits1, attobo.inits2,attobo.inits3, attobo.inits4,attobo.inits5)
set.seed(1128)
attobo.fit <- jags(data = attobo.datjags, inits = attobo.inits, parameters.to.save = attobo.params,
n.chains = 5, n.iter = 100000, n.burnin = 10000, model.file = attobo.model)
attobo.fit.upd <- update(attobo.fit, n.iter =1000)
attobo.fit.upd <- autojags(attobo.fit)
attobo.fit.mcmc <- as.mcmc(attobo.fit)
save(attobo.fit.mcmc, file = "attobo.seekhelp.mcmc.RData")
Convergence Diagnostics
load("attobo.seekhelp.mcmc.RData")
MCMCtrace(attobo.fit.mcmc[,1:4], ind = TRUE, pdf = FALSE )


autocorr.plot(attobo.fit.mcmc[1][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[2][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[3][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[4][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[5][,1:4], ask = FALSE)

attobo.out <- as.data.frame(as.matrix(attobo.fit.mcmc))
## saving data
save(attobo.out, file = "attobo.seekhelp.RData")
load("attobo.seekhelp.RData")
pred.seekhelp <- attobo.out[, grep("lambda[", colnames(attobo.out), fixed = T)]
## estimated count for each participant
pred.seekhelp <- pred.seekhelp[, c(mixedsort(names(pred.seekhelp)))]
strategy.id <- as.numeric(unlist(which(attobo$attribution == "Strategy")))
effort.id <- as.numeric(unlist(which(attobo$attribution == "Effort")))
aptitude.id <- as.numeric(unlist(which(attobo$attribution == "Aptitude")))
pred.seekhelp$strategy.sim <- apply(pred.seekhelp[, strategy.id], 1, mean)
pred.seekhelp$effort.sim <- apply(pred.seekhelp[, effort.id], 1, mean)
pred.seekhelp$aptitude.sim <- apply(pred.seekhelp[, aptitude.id], 1, mean)
mean(pred.seekhelp$strategy.sim)
## [1] 0.2078543
quantile(pred.seekhelp$strategy.sim, probs = c(.025, .5, .975))
## 2.5% 50% 97.5%
## 0.1151397 0.2016900 0.3276598
mean(pred.seekhelp$effort.sim)
## [1] 0.3844671
quantile(pred.seekhelp$effort.sim, probs = c(.025, .5, .975))
## 2.5% 50% 97.5%
## 0.2542277 0.3784910 0.5406089
mean(pred.seekhelp$aptitude.sim)
## [1] 0.7167676
quantile(pred.seekhelp$aptitude.sim, probs = c(.025, .5, .975))
## 2.5% 50% 97.5%
## 0.5287659 0.7127448 0.9281304
error.seekhelp <- attobo.out[, grep("error[", colnames(attobo.out), fixed = T)]
error <- as.vector(unlist(lapply(error.seekhelp, mean)))
mean(error)
## [1] 0.001484177

Expected Counts for each attribution
attobo.posterior.coefs <- select(pred.seekhelp, `strategy.sim`, `effort.sim`,`aptitude.sim`)
names(attobo.posterior.coefs) <- c("Strategy", "Effort", "Aptitude")
attobo.posterior.coefs.long <- gather(attobo.posterior.coefs)
#head(attobo.posterior.coefs.long)
library("ggridges")
ggplot(data = attobo.posterior.coefs.long, aes(x = value, y = key)) +
stat_density_ridges(quantile_lines = TRUE,
quantiles = c(0.025, 0.5, 0.975),
alpha = 0.7)

Bayesian Model - Total Advice
strategy <- ifelse(attobo$attribution == "Strategy", 1,0)
effort <- ifelse(attobo$attribution == "Effort", 1, 0)
total <- attobo$total
attobo.datjags <- list(strategy = strategy, effort = effort, total = total, N = nrow(attobo))
attobo.datjags
Model Specification
attobo.model <- function() {
for(i in 1:N){
# lambda is mean and variance of the poisson distribution
total[i] ~ dpois(lambda[i])
log(lambda[i]) <- b[1] + b[2] * effort[i] + b[3] * strategy[i]
}
#priors
b[1] ~ dnorm(0, 0.0025)
b[2] ~ dnorm(0, 0.0025)
b[3] ~ dnorm(0, 0.0025)
for (i in 1:N){
error[i] <- lambda[i] - total[i]
}
}
attobo.params <- c("b", "lambda", "error")
b.inits <- total.freq$coefficients
attobo.inits1 <- list( "b" = b.inits)
attobo.inits2 <- list( "b" = rep(0, 3) )
attobo.inits3 <- list( "b" = rep(1, 3) )
attobo.inits4 <- list( "b" = rep(-1, 3) )
attobo.inits5 <- list( "b" = rep(2, 3) )
attobo.inits <- list(attobo.inits1, attobo.inits2,attobo.inits3, attobo.inits4,attobo.inits5)
set.seed(1128)
attobo.fit <- jags(data = attobo.datjags, inits = attobo.inits, parameters.to.save = attobo.params,
n.chains = 5, n.iter = 100000, n.burnin = 10000, model.file = attobo.model)
attobo.fit.upd <- update(attobo.fit, n.iter =1000)
attobo.fit.upd <- autojags(attobo.fit)
attobo.fit.mcmc <- as.mcmc(attobo.fit)
save(attobo.fit.mcmc, file = "attobo.total.mcmc.RData")
Convergence Diagnostics
load("attobo.total.mcmc.RData")
MCMCtrace(attobo.fit.mcmc[,1:4], ind = TRUE, pdf = FALSE )


autocorr.plot(attobo.fit.mcmc[1][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[2][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[3][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[4][,1:4], ask = FALSE)

autocorr.plot(attobo.fit.mcmc[5][,1:4], ask = FALSE)

attobo.out <- as.data.frame(as.matrix(attobo.fit.mcmc))
## saving data
save(attobo.out, file = "attobo.total.RData")
load("attobo.total.RData")
pred.total <- attobo.out[, grep("lambda[", colnames(attobo.out), fixed = T)]
## estimated count for each participant
pred.total <- pred.total[, c(mixedsort(names(pred.total)))]
strategy.id <- as.numeric(unlist(which(attobo$attribution == "Strategy")))
effort.id <- as.numeric(unlist(which(attobo$attribution == "Effort")))
aptitude.id <- as.numeric(unlist(which(attobo$attribution == "Aptitude")))
pred.total$strategy.sim <- apply(pred.total[, strategy.id], 1, mean)
pred.total$effort.sim <- apply(pred.total[, effort.id], 1, mean)
pred.total$aptitude.sim <- apply(pred.total[, aptitude.id], 1, mean)
quantile(pred.total$strategy.sim, probs = c(.025, .5, .975))
## 2.5% 50% 97.5%
## 1.642788 1.950798 2.307109
quantile(pred.total$effort.sim, probs = c(.025, .5, .975))
## 2.5% 50% 97.5%
## 2.109232 2.462346 2.855788
quantile(pred.total$aptitude.sim, probs = c(.025, .5, .975))
## 2.5% 50% 97.5%
## 1.857219 2.189896 2.550326
pred.total$diff.1 <- pred.total$strategy.sim - pred.total$aptitude.sim
pred.total$diff.2 <- pred.total$strategy.sim - pred.total$effort.sim
pred.total$diff.3 <- pred.total$effort.sim - pred.total$aptitude.sim
quantile(pred.total$diff.1, probs = c(.025, .5, .975))
## 2.5% 50% 97.5%
## -0.7115193 -0.2327644 0.2565646
quantile(pred.total$diff.2, probs = c(.025, .5, .975))
## 2.5% 50% 97.5%
## -1.008393971 -0.510975172 -0.002368577
quantile(pred.total$diff.3, probs = c(.025, .5, .975))
## 2.5% 50% 97.5%
## -0.2380388 0.2741762 0.8017476
error.total <- attobo.out[, grep("error[", colnames(attobo.out), fixed = T)]
error <- as.vector(unlist(lapply(error.total, mean)))
mean(error)
## [1] 0.001559535

Expected Counts for each attribution
attobo.posterior.coefs <- select(pred.total, `strategy.sim`, `effort.sim`,`aptitude.sim`)
names(attobo.posterior.coefs) <- c("Strategy", "Effort", "Aptitude")
attobo.posterior.coefs.long <- gather(attobo.posterior.coefs)
#head(attobo.posterior.coefs.long)
library("ggridges")
ggplot(data = attobo.posterior.coefs.long, aes(x = value, y = key)) +
stat_density_ridges(quantile_lines = TRUE,
quantiles = c(0.025, 0.5, 0.975),
alpha = 0.7)
